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Abstract 

The Immersed Boundary method is a simple, efficient, and robust numerical scheme for solving PDE in 
general domains, yet it only achieves first-order spatial accuracy near embedded boundaries. In this paper, 
we introduce a new high-order numerical method which we call the Immersed Boundary Smooth Extension 
(IBSE) method. The IBSE method achieves high-order accuracy by smoothly extending the unknown 
solution of the PDE from a given smooth domain to a larger computational domain, enabling the use of 
simple Cartesian-grid discretizations (e.g. Fourier spectral methods). The method preserves much of the 
flexibility and robustness of the original IB method. In particular, it requires minimal geometric information 
to describe the boundary and relies only on convolution with regularized delta-functions to communicate 
information between the computational grid and the boundary. We present a fast algorithm for solving 
elliptic equations, which forms the basis for simple, high-order implicit-time methods for parabolic PDE 
and implicit-explicit methods for related nonlinear PDE. We apply the IBSE method to solve the Poisson, 
heat, Burgers’, and Fitzhugh-Nagumo equations, and demonstrate fourth-order pointwise convergence for 
Dirichlet problems and third-order pointwise convergence for Neumann problems. 

Keywords: Embedded boundary, Immersed Boundary, Fourier spectral method, Complex geometry, 
Partial Differential Equations, High-order 


1. Introduction 

The Immersed Boundary (IB) method was originally developed for the study of moving, deformable 
structures immersed in a fluid, and it has been widely applied to such problems since its introduction [1-3]. 
Recently, the method has been adapted to more general fluid-structure problems, including the motion of 
rigid bodies immersed in a fluid [4] and fluid flow through a domain with either stationary boundaries or 
boundaries with prescribed motion [5, 6]. In this broadened context, we use the term Immersed Boundary 
method to refer only to methods in which (i) the boundary is treated as a Lagrangian structure embedded in 
a geometrically simple domain, (n) the background PDE (e.g. the Navier-Stokes equations) are solved on a 
Cartesian grid everywhere in that domain, and (in) all communication between the Lagrangian structure and 
the underlying PDE is mediated only by convolutions with regularized d-functions. These methods have 
many desirable properties: they make use of robust and efficient Cartesian-grid methods for solving the 
underlying PDE, are flexible to a wide range of problems, and are simple to implement, requiring minimal 
geometric information and processing to describe the boundary. 

The IB method belongs to the broad category of methods known as embedded boundary (EB) methods, 
including the Immersed Interface [7], Ghost Fluid [8], and Volume Penalty methods [9]. These methods 
share a common feature: they enable solutions to PDE on nontrivial domains to be computed using efficient 
and robust structured-grid discretizations; yet these methods differ largely in how boundary conditions are 
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enforced and whether or not the solution is produced in the entirety of a simple domain. Methods which 
compute the solution everywhere in a d-dimensional rectangle admit the simplest discretizations and enable 
the use of high-order discretizations such as Fourier spectral methods. Unfortunately, this simplicity comes 
coupled with a fundamental difficulty: the analytic solution to these problems is rarely globally smooth 
on the entire domain. Consider the one-dimensional Poisson problem A u = f on the periodic interval 
T = [0, 27 t] with Dirichlet boundary conditions u(a) = u(b) = 0 for a b £ T. Even if / £ C°°(T), the 
solution u will typically display jumps in its derivative at the values x = a and x = b (see Figure 2a). 
The lack of regularity in the analytic problem leads to low-order convergence in many numerical schemes, 
including the Immersed Boundary method. 

The advantages of EB methods are substantial enough that significant effort has been expended on 
improving their accuracy [10-26]. Two different approaches are generally taken. The first approach involves 
locally altering the discretization of the PDE in the vicinity of the boundary to accommodate the lack of 
smoothness in the solution. One example of this approach is the Immersed Interface method [7]. Such 
approaches are particularly useful for interface problems where the solution is required on both sides of the 
embedded boundary. When the solution is only needed on one side of the interface, a second approach may 
be taken in which variables are redefined outside of the domain of interest to obtain higher global regularity. 
Improved convergence rates are achieved as a natural consequence of the properties of the discretization 
scheme when applied to smooth problems. Variations of this basic idea have been used by the Fourier 
Continuation (FC) [23-25] and the Active Penalty (AP) [26] methods to provide high order solutions to 
PDE on general domains. 

In this paper, we introduce a new method termed the Immersed Boundary Smooth Extension (IBSE) 
method. This method is a first step in remedying two deficiencies of Immersed Boundary methods: 

1. For generic problems, the IB method produces discretizations that are only first-order accurate in the 
vicinity of domain boundaries (or fluid-structure interfaces) [27]. 

2. The IB method is only able to specify Dirichlet boundary conditions (no-slip, for fluid problems) 
without specialized interpolation schemes subordinate to the geometry [28]. Although this is not of 
interest when solving traditional IB fluid-interface problems, it allows the IB method to be generalized 
for solving other PDE (i.e. reaction-diffusion equations). 

In this paper we restrict our attention to PDE set on stationary domains, and consider only PDE without a 
global divergence constraint. Direct-forcing IB methods produce solutions to PDE that are C°, with jumps in 
the normal-derivative of the solution across the boundary, and converge at first order in the grid-spacing Ax 
[4, 5]. Drawing on ideas from the AP and FC methods, we use Fourier spectral methods to obtain a highly 
accurate discretization, while adding a volumetric forcing to non-physical portions of the computational 
domain to force the solution to be globally C k . We will use the shorthand IBSE-/c to refer to our method 
when we need to explicitly denote the global regularity of the solution that is enforced. High-order accuracy 
is achieved naturally, as a simple consequence of the convergence properties of the Fourier transform. 

The IBSE method retains the essential robustness and simplicity of the original IB method. All com¬ 
munication between the Lagrangian boundary and the underlying Cartesian grid is achieved by convolution 
with regularized (5-functions or normal derivatives of those (5-functions. This allows an absolute minimum 
of geometric information to be used. In the traditional IB method, only the position of the Lagrangian 
structure must be known; the IBSE method additionally requires normals to that structure and an indica¬ 
tor variable denoting whether Cartesian grid points lie inside or outside of the physical domain where the 
PDE is defined. Additionally, since normal derivatives can be accurately approximated by convolution with 
derivatives of regularized (5-functions, Neumann and Robin boundary conditions can be imposed in the same 
way that Dirichlet boundary conditions are in direct-forcing IB approaches. 

In contrast to other approaches based on the smooth extension of the forcing function or the solution 
[19, 23-26], we do not extend the forcing function or solution from a previous timestep. Instead, we smoothly 
extend the unknown solution to the entire computational domain. This approach directly enforces smooth¬ 
ness of the solution, and it allows for high-order implicit- time discretizations for parabolic equations and 
implicit-explicit discretizations for many nonlinear PDE. Remarkably, the coupled problem for the solution 
to an elliptic PDE and that solution’s smooth extension can be reduced to the solution of a relatively small 
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dense system of equations (with size a small multiple of the number of points used to discretize the bound¬ 
ary), along with several FFTs. The dense system depends only on the boundary and the discretization, and 
so it can be formed and prefactored to allow for efficient time-stepping of parabolic equations. 

This paper is organized as follows. In Section 2, we introduce the method, ignoring the details of the 
numerical implementation. The fundamental contributions of this paper are contained in the modification 
to the IB discretization that enforces smoothness of the solution, and in the system of equations (given 
in Section 2.2) that allows for the simultaneous solution of the PDE along with the smooth extension of 
that unknown solution. In Section 3, we discuss the particular numerical implementation that we choose, 
detailing a fast algorithm for solving elliptic equations. In Section 4, we compute solutions to the Poisson 
problem in one and two dimensions, demonstrating high-order pointwise convergence: up to fourth-order for 
Dirichlet problems and third-order for Neumann problems. In Section 5, we discuss discretization of the heat 
equation, detailing a fast algorithm for implicit timestepping, and demonstrate fourth-order convergence in 
space and time. These numerical tests include direct comparisons to the Fourier Continuation [24] and 
Active Penalty [26] methods. Finally, in Section 6, we solve two nonlinear problems: a 2D Burgers’ equation 
with homogeneous Dirichlet boundary conditions, and the Fitzhugh-Nagumo equations with homogeneous 
Neumann boundary conditions. We demonstrate fourth-order convergence in space and time for Burgers’ 
equation and third-order convergence for the Fitzhugh-Nagumo equations. 

2. Methods 

Let SI C C, T = d S4, and E = C \ ft. We will assume that fl is smooth and does not self-intersect. Two 
typical domains are shown in Figure 1. We refer to SI as the physical domain , E as the extension domain , and 
C as the computational domain. We first consider the Poisson problem with Dirichlet boundary conditions: 


A u = f in D, (la) 

u = g on r. (lb) 

For now, we assume that / and g are smooth (C* 00 ) functions defined in f l and F respectively. Let f e be 
a smooth extension of / to C, that is, f e £ C°°(C) and f e (x) = /( x) for all x £ fb The direct forcing 
formulation of the Immersed Boundary method provides a way to solve Equation (1): 


A u(x) + j G(s)5(x — s) ds = f e (x) 
/ u{x)8{x — s) dx = g(s) 

Jc 


in C, 
in r. 


(2a) 

(2b) 


Equation (la) is solved in the entire computational domain C, while the boundary condition is enforced 
by the addition of a singular force G supported on the boundary that acts as a Lagrange multiplier. This 
singular force term leads to jumps in the normal derivative of the solution u across the boundary F. This lack 
of smoothness restricts Immersed Boundary formulations such as Equation (2) to first-order convergence in 
the vicinity of the boundary unless one-sided discrete d-functions are used, which make the implementation 
of the method more difficult [10, 27]. 

We define an example one-dimensional problem to illustrate the limited regularity of solutions produced 
by Equation (2). Let the computational domain be C = T, where the one-dimensional torus T is identified 
with the periodic interval [0,27r], and let the physical domain be given by the interval S2 = T\ [3,4]. The 
extension domain for this problem is E = [3,4]. On this domain, we will solve the problem: 


Au = sin x in SI = T \ [3,4], (3a) 

n = 0 on F = {3,4}. (3b) 


As / is chosen to be / = sin.T in SI, we may choose f e = sin a; for all x £ T. The analytic solution to this 
problem is shown in Figure 2a, along with the solution to the extended problem given by Equation (2). Note 
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(b) Flow around cylinder 


Figure 1: Two different and typical domains to solve PDE on. The physical domain Q is shown in white, 
the extension domain E is shown in gray. In both cases, we have taken the computational domain C to 
be the 2-torus, T 2 = [0, 27t] x [0,27t]. Figure la, shows an irregular, connected domain £1 embedded in C. 
Figure lb shows a domain that would be used to compute flow around a cylinder; here the domain fl is 
itself periodic in one dimension, while the extension region E is not connected. 


that the solution u is only continuous despite the fact that f e £ C°°(T). Choosing f e to be smooth leads 
to low regularity in the solution when the solution to the associated homogeneous problem is nontrivial, as 
is the case here. Without modification, direct discretizations of this problem will exhibit slow convergence 
due to the limited regularity of u in the vicinity of the boundary. 

Rather than choose a smooth extension to the forcing function /, we choose an extension of the forcing 
function that gives a smooth solution on the entire computational domain. Let u e be any smooth extension 
of the unknown solution u into C. We can compute a forcing function F e associated with u e : 

F e = Au e . (4) 

The extended forcing function f e is then defined to be 

fe=Xnf + XeF c , (5) 

where \x denotes the characteristic function of the domain X. Figure 2b shows the extended forcing 
function / e , along with the associated smooth solution u to Equation (3). Since the solution u is smooth, 
we expect that standard discretizations of Equation (2) will converge more rapidly when computed using 
the extension f e than when using the extension f e . We emphasize that the extended forcing function f e 
depends on the unknown solution u. 

This motivates us to define a reformulated version of Equation (2), with the extended forcing function 
f e given by f e as defined in Equation (5): 


A u{x) - (; XEA£ k u)(x ) + J G(s)S(x - s)ds = Xnf(x) 

in C, 

(6a) 

/ u(x)5(x — s) dx = g(s) 

in r. 

(6b) 


c 


Here £ k is any smooth extension operator that satisfies 

£ k : C°(C) n C k ( Q) C k (C), (7a) 

(£ku)(x) = u(x ) \/x £ £l. (7b) 

When restricted to the physical domain H, Equation (6a) reduces to A u = /, the problem of physical 
interest. When restricted to the extension domain E , Equation (6a) reduces to Au = A£ k u. As long as 
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(a) with extended forcing f e 


(b) with extended forcing f e 


Figure 2: The solid lines show the analytic solution u to Equation (3) in the physical domain Cl. The 
dashed line gives the solution u to the extended problem in Equation (2), computed using the extended 
forcing function f e in the left figure, and f e in the right figure. The extended forcing functions are shown 
as the dotted lines. Note the different scales on the forcing functions in the two figures; the forcing is equal 
everywhere in the physical domain Cl. 


u £ C°(C), then Equation (7b) ensures that tt|r = (£ku) |r, and hence u = £kU in E. Again by Equation (7b), 
u = £kU in Cl, and so u = £y-u in C. Since £^u £ C k (C) by Equation (7a), then u £ C k {C). This additional 
regularity of the solution u allows standard discretizations of Equation (6) to converge at a faster rate 
than discretizations of Equation (2). Using C = T d and a simple Fourier-spectral discretization, numerical 
solutions to Equation (6) should converge at 0(Ax k+1 ) [29]. 

In the remainder of this section, we lay out the remaining components needed to fully specify the IBSE 
method that do not depend on the particular choice of numerical implementation. In Section 2.1, we discuss 
how to smoothly extend a known function from Cl to C. In Section 2.2, we outline a system of equations 
that allows us to solve for u and its smooth extension simultaneously. A numerical implementation of this 
method will be discussed in Section 3. 

2.1. Smooth extension of a known function from ft to C 

One way to construct a smooth extension to a function is to solve a high-order PDE. Compared to 
the localized and effectively one-dimensional extension strategies taken in [24, 26], the decision to extend a 
function by solving a fully d-dimensional PDE may appear to be needlessly complex. However, we will show 
in Sections 2.2 and 3 that this choice leads to a robust method that is straightforward to implement and 
requires relatively low computational cost. 

In particular, to compute a C k (C ) extension to a function v £ C k (Cl), we solve the following equation 
in the extension domain E: 

in E, (8a) 

on r, VO < j < k. (8b) 

Here d 3 t;/dn 3 denotes the j tn normal derivative of £ on the boundary T; a computational formula for 
d 3 f/dn 3 is given later in Equation (22). A simple choice of the operator l~L k is the polyharmonic operator 
H k = A fc+1 . From an analytic perspective, this choice is sufficient, although we will show in Section 3.3 
that other choices of H k will produce better results due to numerical issues. The smooth extension £ = £kV 
may then be constructed as 

C = xnv + Xe£,- (9) 

Remark 1. Our choice of smoothed forcing f e only depends on £uu in the extension region E (see Equa¬ 
tion (5)), so there is no need to compute the actual extension £ since the values of £kU in f1 are irrelevant. 


H k £ = 0 
d 3 f d 3 v 
dn 3 ~ dn 3 
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For our purposes, we can think of an ‘extension’ as simply a C k function in the computational domain C 
that shares its first k normal derivatives with w on the boundary T, i.e. the function £. 


As with Equation (1), we may solve Equation (8) in U by extending the problem to all of C. This can 
be accomplished by adding singular forces supported along the boundary: 


^(*) + E(-i) 


k 

E( 

j =0 


. JW dn 3 


(-i y [ 

Jc 


d 3 S(x — s) d 3 v 

- —:- ax = —— r(s) 

dn 3 dn 3 


Vz € T d , 

(10a) 

Vs s r, o < j < k. 

(10b) 


The integral on the left-hand side of Equation (10b) is notation for the action of the distribution d 3 S/dn 3 
on the smooth function £. The boundary condition given in Equation (10b) forces the solution £ to be C k , 
so there is no need to alter this formulation to provide additional regularity. 

For notational convenience, we define the operators T\ and by: 



T *‘« s ) = (/c«»)«(*-<)* 


f c aa 


dS k (x-s) ) T 
dn k ) 


(lla) 

(lib) 


In the language of the Immersed Boundary method, T£ is an interpolation operator; T is a spread operator. 
The notation is suggestive of the fact that these operators obey the adjoint property 

(u,TkF) c = (T k u,F) r (12) 

for all sufficiently smooth u and F, where the notation (•, -) x denotes the L 2 inner product on X. Using 
the operators Tk and T£, we may represent Equation (10) as 


(S r ‘) (« 



(13) 


2.2. A coupled system of equations for u and its extension 

Let the spread (S) and interpolation ( S*) operators be defined as 

SG(x) = J G(s)S(x — s) ds, (14a) 

s*a s ) = f a x )a x ~ s )dx, vssr. (14b) 

Jc 

Notice that S = T 0 and S* =Tq. We can now represent Equation (6) simply as 

Aw - XeA£ + SG = xnf, (15a) 

S*u = g, (15b) 

where £ is defined by the extension equation 

H k f = 0, (16) 


along with the constraint that 
nt = T£u. 
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Equations (15b), (16) and (17) can be combined into one system of equations: 


/ A 

-XeA 

S ) 

n\ 


fxnf\ 

H k T k 


€ 


0 

rp* 

1 k 

rj-i* 

1 k 


F 


0 

V s* 


/ 

\g) 


l 9 / 


(18) 


Equation (18) is the system of equations that allows the IBSE method to solve for u and its extension £ 
simultaneously. The remainder of this paper will be concerned with the discretization and inversion of Equa¬ 
tion (18) and numerical studies demonstrating the accuracy of solutions produced by those discretizations. 
We remark that Equation (18) is equally valid for operators of the Helmholtz type (A — /cl) that appear in 
the discretization of parabolic equations, which will be discussed in Section 5. 


3. Numerical implementation 

For concreteness, we will discuss the numerical implementation in the context of Fourier spectral methods, 
with the computational domain C given by C = T d = [0, 27r] d . We will discretize the domain using a regular 
Cartesian mesh with n points x n discretizing each dimension: Ax = and x n = nAx. Differential 

operators are discretized in the usual way. 

Remark 2. Few of the details depend upon the choice to use Fourier spectral methods, and they are 
chosen because of their simple implementation, computational speed, and high order of accuracy. Inversion 
of the elliptic operator C, as well as the extension operator T-L k is also greatly simplified when using Fourier 
spectral methods. However, any discretization based on a regular Cartesian mesh could be used with minimal 
modifications to the method. 

In order to fully describe a discretization and solution strategy to Equation (18), we must describe several 
key elements. 

1. In Section 3.1, we discretize the spread (5, T k ) and interpolation (S*, T£) operators. 

2. In Section 3.2 we define a new regularization of the 5-function that is accurate and smooth. 

3. In Section 3.3, we define our extension operator H k and show how this choice of W k controls the 
numerical conditioning of the extension problem. 

4. In Section 3.4, we describe an efficient inversion strategy for the IBSE-fc system given by Equation (18). 

5. Finally, in Section 3.5, we provide implementational details and briefly discuss the computational 
complexity of our inversion scheme. 

3.1. Discretization of S, S*, T k , and T k 

Let the boundary T be parametrized by the function X !(s). In all examples in this manuscript, the 
boundary T is one-dinrensional and closed; the single parameter s is defined on the periodic interval [0,27r]. 
The spread operator S : T —> C is defined as 

(SF)(x) = j F(s)S(x-X(s))ds. (19) 

The discrete version of this operator requires a regularized 5-function and a discretization of the integral 
over r. Construction of a regularized 5-function with the properties necessary for the IBSE-fc method is non¬ 
trivial; we will delay discussion of this choice until Section 3.2 and simply denote the regularized 5-function 
as 5. Multivariate 5-functions are computed as Cartesian products of the univariate 5 and also denoted by 
5. Discretization of the integral over T is made by choosing 7ibdy quadrature nodes F = equally 

spaced in the parameter interval [0,27 t] so that X,; = X'(si) and Si = (i — l)27r/nbdy- Quadrature weights 
are computed at each quadrature node to be Ast = || ^-(s»)|| 9 ; this is a spectrally accurate quadrature rule 
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for the integral of smooth periodic functions on T. The discrete spread operator S maps functions sampled 
at points in T to C by 

^■bdy 

(SF)(x) = F{si)S(x - Xi)Asi. (20) 

i=l 

We do not adopt explicit notation to distinguish between the analytic and discretized operators. The number 
of points in the quadrature is chosen so that As ps 2Ax. This choice of node-spacing is wider than that 
recommended for the traditional IB method [1] but has been observed empirically to be the optimal choice in 
other studies of direct-forcing IB methods [4]. The interpolation operator S* may be defined by the adjoint 
property (u,SF) c = (S*u,F) r , but we note that the discrete interpolation operator S* produces a discrete 
function 


(S*u)(sk) = [ u{x)S(x — X(sk)) dx. 

Jc 


( 21 ) 


Discrete integrals over C are straightforward sums computed over the underlying uniform Cartesian mesh. 
Normal derivatives of 6 are computed by the formula [30] 


&>5 

dvJ 


&6 


"■ dr,. 


( 22 ) 


where the Einstein summation convention has been used to indicate sums over repeated indices and d 3 6/dxi 1 ■ ■ ■ 
is computed as Cartesian products of the appropriate derivatives of the one-dimensional 5. For example, 
dS/dn in two dimensions is computed as 


06 ~ ~ 

— = n x o <g) 6 + n v o 8 o . 
on 

Analogous to the definition of S, we define the spread operator for the j th normal derivative, Tjyj, as 

" bdy 


(T(j ) F)(x) = (-1 y 55 nsi)g^(x ^ X,)A S 


(23) 


(24) 


and again define the interpolation operator T U) by the adjoint property (u,T^F) c = (T?^u, i r ) r . Analo¬ 
gous to the action of S *, the operator T U) produces approximations of the j th normal derivative of a function 
defined on C at the quadrature nodes of T. Finally, the composite operator T)* is defined by: 

. T 


n = iT{ 0) t ( * 1} 


[ (fc) 


(25) 


while T k = i T {j)- 


3.2. Construction of a smooth discretization of 5 

The IBSE-A: method is capable of producing solutions that converge at 0{Ax k+1 ). In order to achieve 
this accuracy, our choice of regularized (5-function must satisfy several conditions. For a general value of k: 

1. Its interpolation accuracy must be at least 0(Ax k+1 ). 

2. It must be at least C k , so that its fc th normal derivative is continuous, allowing it to be used in the 
discretization of Tk and T£. 



In addition, we will require that the 5-function has compact (and small) support so that the spread and 
interpolation operators S, S*, Tfc, and T£ may be rapidly applied. The commonly used ‘four point’ 5- 
function is C 1 , has a support of four gridpoints, and produces approximations of S, S *, Tf and T* that are 
accurate to 0( Ax 2 ) [31]. The use of this 5-function regularization is sufficient for the implementation of 
IBSE-1, but it is not sufficient for higher order methods, due to both its limited accuracy and regularity. 

In this manuscript, we discretize the IBSE-/c method for k = 1, 2, and 3, corresponding to second, third, 
and fourth order accuracy in Ax for Dirichlet problems. We therefore need a regularization of 5 that is 
C 3 and has an interpolation accuracy of 0{ Ax 4 ) for smooth functions. We are not aware of any functions 
with these properties currently defined in the literature, and so we construct such a function here. For 
simplicity, we do not attempt to impose other conditions that are often imposed on 5-functions, such as the 
even-odd, condition or sum of squares condition [32]. The strategy that we follow is to choose a 5-function 
with sufficient accuracy but limited regularity and convolve it against itself to increase its smoothness. A 
similar approach was used in [33] to generate a smoother (C 2 ) version of the traditional Peskin four point 
5-function [1] . In order to provide an analytic formula, a base 5-function with a simple functional form must 
be used. We start with the function 

ri-i|r|-|r| 2 + !|r| 3 0 < |r| < 1, 

5/b 4 M = S 1 - ^ M + |r| 2 - | |r| 3 1 < |r| < 2, (26) 

[o 2 < |r|, 

which has 0(Ax 4 ) interpolation accuracy, C° regularity, and a support of four gridpoints [34, 35]. Define: 

5 = 6ib 4 * 5/b 4 * 5ib 4 * 5/b 4 , (27) 

where * denotes convolution. It is clear that 5 £ C 3 , its support is 16 gridpoints, and it is a simple exercise 
to show that convolution preserves interpolation accuracy. The formula for 5 is given in Appendix A. 


Remark 3. While our choice of 5-function regularization limits the spatial discretization in this paper to 
fourth-order accuracy, this is not a fundamental limitation of the method. One way to generalize to higher 
orders would be to construct smoother and more accurate regularizations of the 5-function similar to the 
one that we construct, which maintain finite support for computational efficiency. An alternative approach 
would be to use globally supported regularizations of the 5-function that are C°° (i.e. sine or Gaussian 
functions) and make use of fast sine transforms or fast Gaussian transforms in order to rapidly apply the 
spread and interpolation operators [36, 37]. 

Remark 4. Our construction of 5 is likely not optimal, in that there probably exist C 3 5-function regular¬ 
izations accurate to 0(Ax 4 ) with a support of less than 16 gridpoints (e.g. a 5-function with C 3 regularity 
and accuracy of 0(Ax 2 ) with support of only six gridpoints is defined in [38]). The construction of such a 5- 
function would be worthwhile, allowing for coarser discretizations of problems with closely spaced boundaries 
and faster application of the operators S, S *, X*,, and T£. 


3.3. Choice of extension operator T~L k 

Due to the nullspace of the periodic polyharmonic operator A fe+1 , we choose the operator T-L k used in 
the extension problem given by Equation (8) to be the invertible operator 

U k = A fe+1 + (-1 ) fc+1 0(fc, n). (28) 


Here 0 is a positive scalar function that depends on the smoothness of the extension (fc) and the number 
of Fourier modes ( n ) used in the discretization of the problem. The function 0 is chosen to mitigate the 
numerical condition number of the operator T-L k : 


k = 1 + 


{n/ 2) 2 ( ?c + 1 ) 


(29) 


0 
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Minimizing the condition number n (by taking 0 to be large) must be balanced against the need to resolve 
the intrinsic length scale L = 0~ 1 / 2 ( fc + 1 ) introduced to the problem. In practice, we find that taking 0 to 
be 

/ /n \2(k+l)\ 

0* = max ( 1, ae J J , (30) 

where e is machine precision (2 -52 for double precision and 2~ 112 for quadruple precision) and a = 0.001 
produces stable and accurate solutions. Numerical results are not particularly sensitive to the choice of a 
(see Figure 3a). This stability is demonstrated later (in Section 4.1) with highly accurate solutions shown for 
the IBSE-1, IBSE-2, and IBSE-3 methods up to n = 2 22 . All numerical results in this paper are produced 
using 0* to construct T~L k . 

For large values of n and k , the value of 0 can substantially impact the overall accuracy of the IBSE-fc 
algorithm. We demonstrate this effect in Figure 3 for n = 2 16 , k = 3, and solutions to Equation (3). In 
Figure 3a, we plot the L°° error in the solution u as a function of 0 across 32 orders of magnitude. The 
minimum error observed is 7.63 x 10~ 12 , at 0 = 10 17 , which is within an order of magnitude of 0* (denoted 
by a vertical grey line). Although there is some sensitivity to 0 near the minimum, it is small: the error 
is bounded by 10~ 10 over 6 orders of magnitude of 0. In Figure 3b, we show a zoom of the extension near 
the domain boundary at x = 3 for the solution produced using 0*, which has an error of 1.12 x 10 -11 . The 
extension decays rapidly to 0, but is well resolved by the discretization: there are 150 gridpoints between 
the boundary at x = 3 and the peak at x « 3.014. In Figure 3c, we show a zoom of the extension produced 
when 0 = 10 30 , which gives a solution with an error of 1.23 x 10 -7 , four orders of magnitude worse than the 
error produced when 0 = 0*. Although taking 0 this large leads to a very well-conditioned operator H k , 
the extension decays very rapidly to 0, and the discretization is no longer able to fully resolve the length 
scales: there are only three gridpoints between the boundary at x = 3 and the peak at x ss 3.0004. Finally, 
in Figure 3d, we use 0 = 1. Although there are no fine length-scales to be resolved, ill-conditioning in the 
operator H k leads to an error of 5.2 x 10~ 2 . 


3.f. Inversion of the IBSE-k system (18) 

In this section, we consider solving Equation (18) for an invertible elliptic operator C in place of A. The 
case where C is the periodic Laplacian is complicated by the nullspace of A: even though Equation (18) is 
invertible, the method for inversion that we outline in this section is not directly applicable; details for the 
resolution of this problem are provided in Appendix B. The system of equations is 


/ £ ~xe£ 

s \ 

/ u\ 


(xnf\ 

n k 

T k 

z 


0 

—T£ Tl 


F 


0 

\ s* 

/ 

\c) 


V 9 ) 


This system is large: in two spatial dimensions, the number of equations in the system is 2n 2 + (k + l)?rbdy- 
By block Gaussian elimination, we can find a Schur-complement for this matrix, which we label SC: 


SC = 





(32) 


along with an associated system of equations for the Lagrange multipliers F and G: 



(S*£-' X nf - g\ 

\ T*£-W J' 


(33) 


The size of SC is comparatively small, only (k + l)nbdy square. This equation is the key to the efficiency of 
the algorithm, as it allows the Lagrange multipliers F and G to be computed rapidly without first solving 
for u or £. Once F and G are determined, we may solve for f: 

Z = -{U k )- 1 T kF , 
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(34) 
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(a) L°° error as a function of © 
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Figure 3: The effect that the value 0, used to adjust the conditioning of the extension operator H k , has 
on the solution u of Equation (3). Figure 3a shows a plot of the L°° error against 0. Three special points 
are denoted in this figure: 0 = 0* as defined in Equation (30), denoted by (0); 0 = IO 30 , denoted by (O); 
and 0 = 1, denoted by (D). The solutions corresponding to these special values of 0, along with their 
extensions, are shown in Figures 3b to 3d, respectively. For Figures 3b and 3c, only a zoom of the area near 
x = 3 is shown. See Section 3.3 for a detailed discussion. 
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and once £ is known, u may be computed as 

u = C- 1 (; XD f + Xe£Z ~ SG ). (35) 

The computation of £ and u requires only a series of fast operations: FFTs, multiplies, and applications of 
the spread operators. 

Rapid inversion of the system of equations for F and G given in Equation (33) is not a trivial task. 
Because we have restricted to problems set on stationary domains and the size of SC is small, it is feasible 
to proceed using dense linear algebra. We form SC by repeatedly applying it to basis vectors. This operation 
is expensive: for two dimensions it is 0(N 3 / 2 log N) in the total number of unknowns N = n 2 . The Schur- 
complement is then factored by the LU algorithm provided by LAPACK [39]. Once this factorization is 
computed Equation (33) can be solved rapidly. This Schur-complement depends only on the domain and the 
discretization, so the LU-decomposition can be reused to solve multiple problems on the same domain or in 
each timestep when solving time-dependent problems. See Section 3.5 for a discussion of the computational 
complexity of the method and Table 2 for the actual numerical cost of the method applied to a test problem. 

3.5. Summary of algorithm and numerical complexity 

For a given geometry, choice of elliptic operator £, and discretization size n, the implementation of the 
IBSE-A: method proceeds as follows. 

1. Setup 

(a) The boundary T is discretized as described in Section 3.1, and stencils for the evaluation of 6 and 
its first k normal derivatives are computed at all nodes of the discrete boundary in T to allow for 
the rapid application of the S, S*, T^, and Tf operators. 

(b) The Schur-complement matrix (Equation (32) for invertible £, and Equation (70) for £ = A or 
other non-invertible C) is formed columnwise, by repeatedly applying it to basis vectors. For 
invertible £, this task may be completed as follows (for non-invertible £, the process is nearly 
identical): 

i. One element of either F or G is set to 1, all other elements are set to 0. 

ii. These forces are spread to the Eulerian grid, by the application of T\ and S. 

iii. £, and then u, are computed by applying Equation (34) followed by Equation (35). 

iv. The quantities S*u and —T£u + T££ are computed, and placed in the appropriate rows and 
column of the Schur-complement. 

Note that item (ii) corresponds to application of the right-most matrix in Equation (32), item (iii) 
corresponds to inverting the middle matrix in Equation (32), while item (iv) corresponds to the 
application of the left-most matrix in Equation (32). We remark that although this operation is 
expensive (see Table 2), it is trivially parallelizable, as the task of applying the Schur-complement 
to each individual basis vector is independent. 

(c) The LU decomposition of the Schur-complement is computed. 

(d) Optionally, these items may be saved, bypassing the expensive computation and factorization of 
the Schur-complement in future uses of the IBSE-fc method when applied to the same geometry, 
£, and discretization size n. 

2. Solve 

(a) The right-hand side of Equation (33) (or equivalently, Equation (70) when C is not invertible) is 
computed. 

(b) Using the pre-computed LU decomposition of the Schur complement matrix, the Lagrange mul¬ 
tipliers F and G are found by solving Equation (33) (or equivalently, Equation (70) when C is 
not invertible). 

(c) The smooth extension £, and finally u, are computed via Equations (34) and (35) (or equivalently, 
Equation (72) when £ is not invertible) 
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Dimension 

Cost of FFT 

Cost of Setup 

Cost of Solve 

2 

O(NlnN) 

0(N 3 / 2 In TV) 

O(NlnN) 

3 

O(NlnN) 

0(N 2 ) 

0(N 4 / 3 ) 


Table 1: Numerical cost of the IBSE-fc algorithm as a function of N = n d for dimensions d = 2 and d = 3, 
along with the cost of an FFT. The ‘Cost of Solve’ refers to the cost of solving Equation (18) once, which is 
equivalent to the cost of solving one elliptic problem or taking one implicit timestep of a parabolic problem. 


In items 1(b), 2(a), and 2(c), application of the FFT and IFFT is used to move variables between frequency 
and physical space, while differential operators are applied in frequency space and multiplications with 
characteristic functions are carried out in physical space. 

The computational work for this method can be broken into two main tasks: (i) the setup items, 
listed above, which are dominated by the formation and factorization of the Schur-complement defined in 
Equation (32) and (ii) the production of one solution u given the Schur-complement, its factorization, the 
body forcing /, and the boundary conditions g. We summarize the scaling of the algorithm in Table 1; 
see Table 2 in Section 5.1 for the actual numerical cost of the method applied to a test problem. In two 
dimensions, the cost of a solve scales like the cost of an FFT; in three dimensions the cost is slightly 
worse: 0{N 4 ^ 3 ) in the total number of unknowns N = n 3 . The higher cost of the algorithm in three- 
dimensions is due to the cost of factoring and inverting the Schur-complement, however, these matrices are 
highly structured. It is quite possible that the cost of inversion could be reduced with an appropriately 
preconditioned iterative method, by directly exploiting the structure of the matrix, e.g. with an inverse 
Multipole Method [40], or by indirectly exploiting that structure using an algorithm like HODLR [41] or 
Algebraic Multigrid [42]. 

4. Results: Poisson equation 

4-1- One-dimensional test problem 

To demonstrate both the power and the limitations of the IBSE method, we compute the solution to 
the one-dimensional example defined in Equation (3). Solutions are computed on grids ranging from n = 2 4 
to 2 22 points, in both double- and quadruple-precision arithmetic, for the IB, IBSE-1, IBSE-2, and IBSE-3 
methods. Figure 4 shows a refinement study demonstrating the expected first-, second-, third-, and fourth- 
order accuracy in L°° for the IB, IBSE-1, IBSE-2, and IBSE-3 methods, respectively. The IB method 
achieves an error of 4.23 x 1CF 7 with n = 2 22 grid points. Imposing additional smoothness on the solution 
allows this error to be matched by the IBSE-1 method at n = 4096, by the IBSE-2 method at n = 1024, and 
by the IBSE-3 method at n = 512. For the IBSE-3 method, this is a factor of nearly 10000 less gridpoints; 
equivalent to a factor of 100 million less gridpoints for two-dimensional problems. In practice, obtaining 
solutions accurate to six digits with the traditional IB method is often impractical; with the IBSE method 
this kind of accuracy can be achieved on reasonably sized grids. 

In double precision, all of the methods achieve the expected convergence rate up to some value of n. The 
double-precision solutions begin to diverge from the quadruple-precision solutions at n = 2 1 ' for IBSE-1, at 
n = 2 16 for IBSE-2, and at n = 2 12 for IBSE-3. Using 0 = 0* (see Section 3.3) to control the condition 
number of the extension operator TL k provides sufficient numerical stability at all values of n to enable the 
computation of solutions accurate to 10-11 digits. In quadruple precision, all methods exhibit the expected 
order of convergence across a wide range of values of n, although the fourth order method begins to fall 
short of the expected path of convergence for extremely fine grids {n > 2 20 ). 

Remark 5. We note that this test problem, as well as the test problems in Sections 5.1, 6.1 and 6.2, are 
exterior problems set on periodic domains. Periodicity of solutions in the physical domain f2 is enforced 
naturally in these problems through the use of the Fourier basis. 
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Figure 4: L°° error for solutions to Equation (3) produced with IB (0), IBSE-1 (A), IBSE-2 (O), and IBSE- 
3 (D), demonstrating O(Ax), 0(Ax 2 ), 0(Ax 3 ), and 0( Aa; 4 ) convergence, respectively. Solutions computed 
in double precision (64-bit) are shown with solid lines; solutions computed with quadruple precision (128-bit) 
are shown with dashed lines. 


4-2. Two-dimensional test problem: solution inside a circle 

Let Cl = be the circle of radius 2 centered at (7r,7r), and identify C = T 2 with the square 

[0, 27t] x [0, 2n]. We will solve the Dirichlet problem: 

Art = —4 
u — 0 

which has an analytic solution u a given by 

u a = 4- (x - 7r) 2 - (y - 7r) 2 . (37) 

We solve this problem in double-precision arithmetic with the IBSE-1, IBSE-2, and IBSE-3 methods for 
n = 2 4 to n = 2 11 . In Figure 5, we show L°° error as a function of n, indicating second-, third-, and 
fourth-order accuracy for IBSE-1, IBSE-2, and IBSE-3 respectively. 

Using the solutions to this problem, we demonstrate the property that the IBSE method produces globally 
smooth solutions. In Figure 6, we show the solution u to Equation (36), along with its first, second, and third 
derivatives in the ^-direction, produced using the IB, IBSE-1, IBSE-2, and IBSE-3 methods with n = 2 9 . 
For simplicity, all functions are shown only along the slice y = 7r, and the intersection of this slice with the 
boundary T is shown as gray vertical lines. Derivatives are computed spectrally using the FFT. We expect 
solutions produced by the IBSE-fc method to exhibit discontinuities in the (k -I- l) st normal derivative across 
interfaces, and for all derivatives of lower order to be continuous. We see this expectation validated in the 
solutions presented in Figure 6: functions shown in plots with shading behind them in Figure 6a are at least 
continuous , all others are discontinuous. This property of global smoothness of the solutions enables the 
Fourier-spectral discretization to obtain high-order accuracy. 

4-3. Two-dimensional test problem with Neumann boundary conditions 

Traditional Immersed Boundary approaches are unable to provide solutions to Neumann problems since 
convolution-style estimation of the normal derivative is not convergent at the boundary due to the lack 


in Cl, (36a) 

on F = dCl, (36b) 
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Figure 5: L°° error for solutions to Equation (36) produced with IBSE-1 (A), IBSE-2 (O), and IBSE-3 
(D), demonstrating 0(Ax 2 ), 0(Ax 3 ), and 0(Ax 4 ) convergence, respectively. 



x 


(a) (b) 

Figure 6: A demonstration of the smoothness of the solution u to Equation (36) produced by the IB, IBSE-1, 
IBSE-2, and IBSE-3 methods. Figure 6a shows the solution u and the solution’s first, second, and third x 
derivatives. All plots show only the slice y = tv from x = 0 to x = 27r and vertical gray lines denote the 
intersection of that slice with the boundary T at x = 7r — 2 and x = tt + 2. The functions in plots with 
shading behind them (to the lower left) are at least continuous , all other functions shown are discontinuous. 
Solutions produced using the IBSE-/c method exhibit discontinuities in the (fc + l) st normal derivative across 
T; all lower-order derivatives are continuous. Scales of the three plots in the upper right have been omitted. 
For these plots, zooms of the region near the boundary at x = tt — 2 are shown in Figure 6b. The x-axis in 
these plots runs from x = tt — 2 — 0.3 to x = tt — 2 + 0.3. 
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of regularity of the solution. The IBSE method does not have this limitation because of the additional 
smoothness of the solution. The only modification necessary to adapt Equation (18) to solve Neumann 
problems is to change the S and S* operators that enforce and specify the boundary conditions. Consider 
the general Neumann problem: 


in f2, 
on T. 

Define 

T W F(x) = j/M d5{X d ~ S) ds, 

Tfotiis) = -1 g(s ) d6(x g ~ s) dx, Vs e r. 


A u = f 
du 


dn 


= 9 


(38a) 

(38b) 


(39a) 

(39b) 


The analogue to Equation (18) is 


/ A 

-Xe A 

n k 

T w\ 

u\ 
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( Xnf\ 
0 

rn* 

1 k 

W) 

rn* 

1 k 
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F 

\G) 


0 

\ 9 J 


(40) 


Remark 6. It is just as simple to modify the IBSE method to allow the solution of Robin problems. 
Rather than replacing the upper-right and lower-left operators in Equation (18) with Tjq) and 7)* .. they 
should instead be replaced with linear combinations of S, Tqp S*, and T* ± y 


Let = B\ ((7r, 7r)). We will solve the problem 


A u = e smx (cos“ x — sin a;) — cos y 


du 

dn 


= cos 2 x e sin x — sin 2 y 


in f2, 
on T, 


(41a) 

(41b) 


which has the analytic solution 


u a = e 


cos y. 


(42) 


We solve this problem with the IBSE-1, IBSE-2, and IBSE-3 methods, for which we expect first-, second-, 
and third-order convergence, respectively. The reason for the lower order of convergence than in Dirichlet 
problems is that our discretization of 7)* , is one order less accurate than S* when acting on functions of the 
same regularity. Despite this, the solution u produced by IBSE-fc still has global C k regularity. Solutions 
are computed in double precision, for n = 2 4 to n = 2 11 . Figure 7 shows L°° error as a function of n. 
Second- and third-order convergence is observed for the IBSE-2 and IBSE-3 methods. The IBSE-1 method 
converges at a rate that is slightly less than first-order. 


5. 


Results: heat equation 

We now consider solving the heat equation with Dirichlet boundary conditions: 

Ut — A u = f(x,t) in D, 

u = g(x , t ) on T. 


(43a) 

(43b) 
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Figure 7: L°° error for solutions to Equation (41) produced using IBSE-1 (A), IBSE-2 (O), and IBSE- 
3 (D) demonstrating slightly less than first-order convergence, second-order convergence, and third-order 
convergence, respectively. 


Since the IBSE method can directly invert elliptic equations, it can easily be adapted to provide higlr-order 
implicit-time discretizations for the heat equation. Many implicit methods could be used; we choose BDF4 
[43]: 

(I - ^AtA) u n+1 = (l2At/(t + At) + 48 u n - 36U"- 1 + 16 m ”" 2 - 3 u n ~ 3 ) . (44) 


Applying this scheme to Equation (43b) is equivalent to solving the problem: 

Cu n+1 = / in Cl, (45a) 

u n+1 = g{x,t + At) on T, (45b) 

where we define: 

H = I — — At A, (46a) 

25 

/ = — (12At/(t + At) + 48 u n - 36U”- 1 + 16u"- 2 - 3u”- 3 ) . (46b) 

Zo 


This may be solved using the algorithm described in Section 3.4. There are several comments worth making 
regarding this discretization: 

1. The use of an implicit scheme allows for large timesteps to be taken. For all test problems presented 
in this section, we choose At to be: 


At = 


^final 

~Y~ 



(47) 


where [•] is the ceiling function, so that At « Ax/2. 

2. The fact that C depends on At implies that the Schur-complement depends on At as well. This means 
that the Schur-complement must be reformed and refactored whenever At is changed, complicating 
the use of a method that uses adaptive timestepping. 

3. Modifying this formulation to solve the heat equation with Neumann and Robin boundary conditions 
may be done in the same way as shown for the Poisson equation in Section 4.3. 
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n x n 

IBSE-1 

IBSE-2 

IBSE-3 

timesteps to 
t = 0.1 

prep time 

timestep 

prep time 

timestep 

prep time 

timestep 

32 x 32 

674 

39 

895 

39 

1451 

43 

2 

64 x 64 

957 

34 

1776 

39 

2371 

39 

3 

128 x 128 

1310 

30 

2207 

31 

3241 

37 

5 

256 x 256 

1833 

21 

2833 

24 

3917 

23 

9 

512 x 512 

3478 

21 

4493 

22 

5782 

20 

17 

1024 x 1024 

5328 

19 

7521 

19 

10187 

19 

33 

2048 x 2048 

9951 

18 

13707 

18 

18686 

18 

66 


Table 2: The computational time, normalized by the time to compute an FFT, required for precomputation 
(labeled ‘prep time’), and per timestep of Equation (48) (labeled ‘timestep’). We note that even for n = 2048, 
the ‘prep time’, shown here as requiring 18686 FFTs, requires only 50 minutes of wall time using serial code 
that has not been carefully optimized. See text for further discussion. The number of timesteps taken to 
advance Equation (48) from t = 0 to t = 0.1 is also shown. 


5.1. Solution in a periodic domain outside an obstacle 

In this section we demonstrate high-order convergence to a heat-equation set on a periodic domain with 
a simple obstacle in it. To provide direct numerical comparison with results from the Active Penalty (AP) 
method, the following problem is from Section 6.2 of [26]: 

u t — A u = [cosy + e slnx (sinx — cos 2 x)] cos t — (e slna: + cosy) sint in ft, (48a) 

u(x,y,t. = 0) = e sln x + cos y in fi, (48b) 

u = (e sln x + cos y) cos t on T, (48c) 

where the physical domain is fl = T 2 \ B 1 / i (ir, n). Equation (48) has the analytic solution 

u a = (e sm x + cos y) cos t. (49) 

The initial condition is integrated from t = 0 to t = 0.1 by solving Equation (45b) at each timestep. Startup 
values for BDF4 are computed by evaluating the analytic solution at t = —At, t = —2A t, and t = —3A t. 
In Figure 8, L°° errors at t = 0.1 for solutions generated with the IBSE-1, IBSE-2, and IBSE-3 methods are 
shown, demonstrating second-, third-, and fourth-order convergence in space and time, respectively. On this 
plot, we also show errors from the second- and third-order AP methods with n = 512, the finest resolution 
reported in [26]. Our method yields errors lower by several orders of magnitude. 

In contrast to the AP method [26], we are able to use an implicit time discretization, enabling the 
use of large timesteps. The ability to take large timesteps efficiently comes at the cost of a significant 
precomputation for the IBSE method, a cost that the AP method does not incur. This precomputation 
prevents the IBSE method as described in this paper from being efficiently applied to moving boundary 
problems, a limitation not faced by the AP method. 

In Table 2, we provide the computational time required by the IBSE-1, IBSE-2, and IBSE-3 methods 
to setup the problem (dominated by the formation and factorization of Equation (32)), the time required 
to take one timestep, and the number of timesteps taken to advance the equation from t = 0 to t = 0.1. 
Computational times are given as a multiple of the time required to take one FFT. In our code, 8 FFTs 
are used per timestep regardless of k: the time required to take a timestep does not vary significantly 
across the methods. The precomputation time increases approximately linearly in fc, but even for the finest 
discretization presented, the precomputation time is not prohibitive: for n = 2048 and the IBSE-3 method, 
the precomputation requires about 50 minutes using serial code that has not been carefully optimized. 
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Figure 8: L°° error for solutions to Equation (48) produced using IBSE-1 (A), IBSE-2 (O), and IBSE-3 
(□), demonstrating 0(Ax 2 ), 0(Ax 3 ), and 0( Ax 4 ) convergence, respectively. Results from [26] are shown 
for n = 512, for the second- (A) and third-order(®) Active Penalty methods [26]. 


5.2. Solution inside a parametrically defined region 

In this section, we demonstrate high-order convergence to a heat equation set inside a parametrically 
defined domain. This problem is from Section 6.1 of [24], allowing direct comparison of the IBSE method 
with the Fourier Continuation method [24] . For convenience, we rescale the domain of the problem. We will 
solve 


u t — Au = 7r (2 — A<j>) cos(7 T(j)) + 7r 2 (</> 2 + </>()) sin(7r</>) in fl, (50a) 

u(x, y, t = 0) = sm(TT(j>(x, y,t = 0)) in f 2, (50b) 

u = sin(7r</>) on T, (50c) 

where 

(p(x,y,t) = 9 + lj +4^^-^ + l^ +2 1. (51) 

The domain Q is the region inside the boundary F defined by the parametric equations: 

x{9) = (lOsin 2 (20) + 3cos 3 (20) + 40) cos0/20 + 7r, (52a) 

2/(0) = (lOsin 2 (20) + 3cos 3 (20)+4O)sin0/2O + 7r, (52b) 


for 0 < 0 < 27 t. Equation (50) has the analytic solution u = sin(7n j>). The initial condition is integrated from 
t = 0 to t = 0.01 by solving Equation (45b) at each timestep. Startup values for BDF4 are computed by 
evaluating the analytic solution at t = —At, t = —2At, and t = —3A t. In Figure 9, we report the maximum 
error over space and time, demonstrating second-, third-, and fourth-order convergence in space and time 
with the IBSE-1, IBSE-2, and IBSE-3 methods, respectively. Figure 9b shows a solution to Equation (50), 
together with its C 1 extension, computed using IBSE-1 and n = 2 9 at t = 1.0. Also shown in Figure 9a are 
the results produced by the FC method from [24]. The errors for the FC method are reported using Ax in 
[25]; in Figure 9 we plot the results as a function of n so that the values of Ax are equivalent (accounting 
for the rescaling of the problem). For both the IBSE and FC methods, the timestep has been taken low 
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Figure 9: Figure 9a shows L°° errors for Equation (50) at t = 1.0 produced by the IBSE-1 (A), IBSE-2 
(O), and IBSE-3 (D) methods, demonstrating 0(Ax 2 ), 0(Ax 3 ), and 0(Ax 4 ) convergence, respectively. 
The dashed line shows results from the Fourier Continuation method [24]; the value of ‘n’ at which we plot 
their results is set so that the grid-spacing Aa; is the same for the FC and IBSE methods, taking into account 
our rescaling of the domain. Figure 9b shows the numerical solution to Equation (50) computed with the 
ISBE-1 method and n = 2 9 at t = 1.0. The boundary T = dfl of the domain is shown in black. Note that 
the physical domain f 1 is contained inside the boundary. 


enough so that the dominant error is the spatial error. Despite the fact that the FC method is fifth-order 
in space, the IBSE-3 method is able to achieve comparable errors at the finest grid-spacing reported 1 . For 
all of the FC results, the timestep is taken to be 1 x 10 -5 . Although their Alternating-Direction Implicit 
(ADI) scheme is unconditionally stable, it is only second-order accurate, and thus requires the use of small 
timesteps to obtain errors small enough to not impede the rapid convergence of the spatial discretization. 
The ADI nature of their scheme complicates implementation of higher-order schemes, although this issue 
may be resolved using high-order Richardson extrapolation [23]. The structure of the IBSE method allows 
for the use of simple and high-order timestepping methods. In these computations, we use BDF4 and are 
able to take substantially larger timesteps while still obtaining comparable error. 


6. Results: nonlinear problems 

One significant advantage of the IBSE method is the simplicity with which the elliptic and parabolic 
equation solvers can be integrated into methods to solve more complicated PDE. We solve two nonlinear 
problems, the Burgers’ equation and the Fitzhugh-Nagumo equations, as a demonstration. The IBSE method 
is able to obtain high-order convergence in both space and time using simple timestepping methods and a 
straightforward pseudo-spectral computation of nonlinearities. 


1 The method used in [25] was based on an oversampled formulation near the boundary. Newer FC methods [44] do not 
require this oversampling, effectively halving the number of discretization points in each dimension required to obtain a given 
accuracy. 
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Solutions 

L°° difference 

ratio 

log 2 ratio 

2 6 , 2 7 

3.65e-02 

- 

- 

2 7 , 2 8 

9.75e-03 

3.7 

1.9 

2 s , 2 9 

6.87e-04 

14.2 

3.8 

2 9 , 2 10 

2.83e-05 

24.3 

4.6 


Table 3: Unnormalized L°° differences between solutions to Equation (53) at successive levels of grid re¬ 
finement, computed using IBSE-3. Ratios of the L°° differences are computed, indicating fourth-order 
convergence in space and time. 


6.1. Burgers’ equation 

In this section, we apply the IBSE method to a two-dimensional homogeneous Burgers’ equation, a 
simple model with a nonlinearity that has the same form as Navier-Stokes: 

Ut — vAu + u ■ Vu = 0 in U, (53a) 

u = 0 on E, (53b) 

u{x,y,t = 0) = ip(x,y), (53c) 

Here U is the region outside of the polar region defined by r = 1 + cos(f? + 7r/4)/4, shifted by (37r/2, 37t/2); 
this domain is shown in Figure 10a. The initial condition is given by the Gaussian pulse if(x,y) = 
2 e -40(>-2.5) g—40(y—4.3) ^ anc [ v _ q. 01. As with the other examples, the computational domain C is taken 
to be T 2 = [0,2?r] x [0,2 tt], 

We integrate this equation in time using a fourth-order implicit-explicit (IMEX) Backward Differentiation 
formula [45]: 

^ u n+1 - 4 u n + 3b"" 1 - \u n ~ 2 + -u n ~ 3 = 

12 3 4 

At [.I(u n+1 ) + 4 £(u n ) - 6£(u n ~ 1 ) + 4 £{u n ~ 2 ) - £{u n ~ 3 )\ , (54) 

where X evaluates the stiff terms in ut and £ evaluates the non-stiff terms in ut- For this problem, 
X{ii) = vAu and £{u) = —u ■ \7u. In practice, timestepping involves the explicit computation of a forcing 
function: 

= 4 u n - 3 u n ~ x + | u n ~ 2 - ^ u n ~ 3 + At [4 £{u n ) - 6 £{u n ~ 1 ) + 4 £(u n ~ 2 ) - £{u n ~ 3 )} , (55) 

followed by solution of the equation 

in fl, (56a) 

on T. (56b) 

Since u n is smooth in the entire computational domain C , the explicit terms such as £(u n ) may be accurately 
computed using simple pseudo-spectral methods. Equation (56) is solved using the method described in 
Section 3.4. For this test, the startup values n -1 , u~ 2 , and u~ 3 are all taken to be if. We integrate 
Equation (53) to t = 2 with a timestep of At = Aa:/20, for n = 2 6 to n = 2 10 , using IBSE-3 when solving 
Equation (56). Figure 10a and Figure 10b show the initial condition and a zoom of the solution at t = 2 
for n = 2 10 . Convergence is assessed by comparing the ratio of the L°° difference between solutions at 
different resolutions, computed over the spatial locations in the coarsest grid. Results are shown in Table 3, 
demonstrating fourth-order convergence in both space and time. 


I -^i/AtA )u n+1 =f 


n+1 


= o 


21 






Figure 10: The initial condition, geometry, and solution at t = 2.0 to Equation (53). The extension domain 
E is indicated by the hatched region. The physical domain fi is everything outside of E. 


6.2. Fitzhugh-Nagumo equation 

The Fitzhugh-Nagumo equation [46] is a system of nonlinear reaction-diffusion equations often used to 
model excitable media in biological systems (see e.g. [47]). We solve the Fitzhugh-Nagumo equations with 
homogeneous Neumann boundary conditions and an initially unperturbed recovery variable w: 


vAv + v(a — u)(l — v) + w = 0 

in f2, 

(57a) 

wt + e(yw — v) = 0 

in O. 

(57b) 

dv 



dn. = ° 

on r 

(57c) 

C5 

C-K 

II 

O 

II 

2L 


(57d) 

jf^ 

C"-K 

II 

O 

II 

O 


(57e) 


The physical domain Q and the computational domain C are the same as those used for the solution to 
Burgers’ equation in Section 6.1. The initial data is the pulse ip{x,y) = 2e~ 100 ^ x ~ 2 5 ^ e -10 ^ -4 ' 3 ^ . High- 
order time integration is implemented using the same IMEX-BDF scheme used in Section 6.1, with the 
implicit terms given by 



and the explicit terms given by 

g f v\ _ fv{a — v)(v — 1) — w 

\w J e(i> — 7 w). 

We set the parameters to be a = 0.1, 7 = 2, e = 0.005, and v = 0.001. For this test, the startup values v~ l , 
u~ 2 , and v -3 are all taken to be ip. We integrate Equation (57) to t = 200 with a timestep of At = Ax/2, 
for n = 2 6 to n = 2 10 with the IBSE-3 method when solving for the updated value of v. Figure 11 shows the 
initial condition ip, as well as the solutions v and w at time t = 200. Convergence is assessed by comparing 
the ratio of the L°° difference between solutions at different resolutions, computed over the spatial locations 
in the coarsest grid. Results are shown in Table 4, demonstrating at least third-order convergence in both 
space and time, as expected when solving with Neumann boundary conditions. 
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Figure 11: The initial condition for v , and the solutions to v and w at t = 200 for Equation (57) generated 
with the IBSE-3 method. The extension domain E is indicated by the hatched region. The physical domain 
12 is everything outside of E. 


Solutions 

L°° difference 

ratio 

log 2 ratio 

Solutions 

L°° difference 

ratio 

log 2 ratio 

2 6 , 2 7 

2.23e-01 

- 

- 

2 6 , 2 7 

1.63e-02 

- 

- 

2 7 , 2 8 

2.38e-02 

9.4 

3.2 

2 7 , 2® 

2.25e-03 

7.3 

2.9 

2 s , 2 9 

1.86e-03 

12.8 

3.7 

2 8 , 2 9 

1.91e-04 

11.8 

3.6 

2 9 , 2 10 

2.68e-04 

6.9 

2.8 

2 9 , 2 10 

8.68e-06 

22.0 

4.5 


(a) Refinement study for v (b) Refinement study for w 


Table 4: Unnormalized L°° differences between solutions to Equation (57) at successive levels of grid refine¬ 
ment, computed with IBSE-3. Ratios of the L°° differences are computed, indicating at least third-order 
convergence in space and time. 
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7. Discussion 


We have developed a new numerical method, the IBSE method, for solving elliptic and parabolic PDE 
on general smooth domains to an arbitrarily high order of accuracy using simple Fourier spectral methods. 
Accuracy is obtained by forcing the solution to be globally smooth on the entire computational domain. 
This is accomplished by solving for the smooth extension of the unknown solution, enabling the method to 
directly solve elliptic equations. Remarkably, this computation can be effectively reduced to a small, dense 
system of equations, and the cost to invert this system can be minimized by precomputing and storing its 
LU-factorization. This provides a very efficient algorithm for implicit timestepping on stationary domains. 
The method requires minimal geometric information regarding the boundary: only its position, normals, and 
an indicator variable denoting whether points in the computational grid are inside of the physical domain 
Q or not, enabling simple and robust code. The method is flexible enough to allow high-order discretization 
in space and time for a wide range of nonlinear PDE using straightforward implicit-explicit timestepping 
schemes. 

The IBSE method shares some similarities with Fourier Continuation (FC) [23-25] and Active Penalty 
(AP) [26] methods. The FC method uses Fourier methods to obtain a high-order discretization, and the 
idea that smooth, non-periodic functions can be turned into smooth and periodic functions by extending 
them into a larger domain in some appropriate way. There are two key differences between the FC method 
and the IBSE method. The FC method (i) uses dimensional splitting to reduce the problem to a set 
of one-dimensional problems and (ii) relies on data from the previous timestep in order to generate the 
Fourier continuations that allow for their high-order spatial accuracy. This forces the FC method to use 
an Alternating-Direction Implicit scheme to take large stable timesteps when solving parabolic equations, 
complicating the implementation of high-order timestepping. In addition, the FC method requires the use 
of an iterative solver for computing solutions to the Poisson equation, complicating an efficient discretiza¬ 
tion of the incompressible Navier-Stokes equations, although the FC method has been used to solve the 
compressible Navier-Stokes equations to high order [23]. In contrast to the conditioning issues faced by the 
IBSE method, the one-dimensional nature of the Fourier-continuation problem allows the FC method to 
use high-precision arithmetic in certain precomputations to effectively eliminate the conditioning problem 
inherent in constructing smooth extensions [48] . This enables the FC method to obtain stable methods that 
converge to higher order than can be achieved by the IBSE method in double-precision arithmetic. 

The AP method, like the FC method, relies on data from previous timesteps in the way that it imposes 
smoothness on the solution of the PDE. A large drag force is applied that penalizes deviations of the solution 
in the extension region from the smooth extension of the solution at the previous timestep. This dependence 
on data from previous timesteps forces the AP method to use explicit timestepping when solving parabolic 
equations. The smoothness constraint is also enforced only approximately ; solutions are C 1 regardless of how 
many derivatives are matched when generating the extension functions. We believe this property explains 
both the relatively small difference in L°° error between the second-order AP method and the second-order 
IBSE-1 method, and the orders of magnitude difference in L°° error between the third-order AP method 
and the third-order IBSE-2 method observed in solutions to Equation (48). Despite these disadvantages, 
the AP method does not require the solution of a dense system of equations and should be immediately 
applicable to moving boundary problems. 

We have only implemented two-dimensional examples in this paper. The IBSE method extends to three 
dimensions without any changes, although there is one minor difficulty that must be resolved: the production 
of an accurate enough quadrature for the discretization of the boundary integrals in the spread operators 
S and T*,. For two-dimensional problems, this comes nearly for free since the simple quadrature rule given 
in Section 3.1 is spectrally accurate for closed one-dimensional curves. High-order quadrature rules for two- 
dimensional surfaces are more complicated but well-developed, and high-order surface representation has 
been incorporated into the Immersed Boundary method [49, 50]. 

We have discretized the IBSE-fc method for 1 < k < 3. In this work, we are limited to k = 3 largely by 
the accuracy and regularity of our regularization of S. This limitation is not fundamental. There are at least 
two paths forward to obtain higher-order accuracy. Smoother and more accurate ^-functions with compact 
support could be constructed, however the method that we use to construct S (convolving 5-functions with 
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limited regularity against themselves) would most likely not be of use here. For example, to discretize 
the IBSE-5 method, we could start with the sixth-order analog to the base (5-function that we use (see 
Section 3.2), which has a support of six grid points [35]. To obtain a C 5 5-function, this would need to be 
convolved against itself 5 times, resulting in a 5 with a support of 30 gridpoints (900 for a two-dimensional 
5-function). An alternate approach is to use non-compactly supported 5-functions, combined with fast 
transform methods for the application of the spread and interpolation operators [36, 37]. 

When solving the IBSE-A: equations for high values of k, the conditioning of the extension operator T~L k 
given in Equation (7) becomes problematic. Although the conditioning can be controlled by an appropriate 
choice of 0, the length scales introduced by 0 rapidly become unresolvable as the discretization is refined. A 
different choice of extension operator, e.g. a polynomial of the Laplacian whose spectra is explicitly designed 
to minimize the condition number while limiting the introduction of fine scales, could potentially be used 
to mitigate this problem. 

In this manuscript, we have developed two of the necessary components for the implementation of a 
high-order numerical scheme for the incompressible Navier-Stokes equations based on the IBSE method: 
a direct solver for the Poisson equation and a high-order implicit-explicit discretization of a nonlinear 
advection-diffusion equation. The missing component is how to impose a global divergence constraint without 
compromising the accuracy of the discretization. This is a non-trivial problem even for finite-difference 
discretizations on curvilinear domains [51]. We have obtained preliminary results for solving the Stokes 
equation using the IBSE method; this work will be presented in a forthcoming contribution. 

Finally, the high-efficiency that we obtain is only currently achievable on stationary domains. This is 
because the inversion method that we use in this paper for the IBSE method relies on an expensive pre- 
computation that depends on the physical domain and the discretization. However, there is no fundamental 
obstacle to the application of the IBSE method to problems with moving domains. Instead, the challenge 
is to find a robust and efficient method to invert the IBSE system in Equation (18) that does not require 
substantial precomputation. Recent progress has been made for preconditioning similar systems of equations 
for the simulation of rigid-body motion in an Immersed Boundary framework [4]. The integration of these 
ideas into the IBSE method to allow for simulation of moving boundary problems will be an area of active 
future research. 
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Appendices 

A. Formula for C 3 ^-function with interpolation accuracy of 0( A® 4 ) 


We give the formula here for S , defined in Section 3.2. For brevity, we define <5 for r > 0, which fully 
defines the function as 6 is even. The formula for S is defined piecewise; each portion is a fifteenth order 
polynomial. Coefficients are provided in the table below. 



0 < r < 1 

1 < r < 2 

2 < r < 3 

3 < r < 4 

4 < r < 5 

5 < r < 6 

6 < r < 7 

7 < r < 8 

r >= 8 

T° 

12949745023 

3177441629 

21914742667 

4094824493 

-1606651889 

163201885541 

1005507698627 

-5005877248 

0 


20432412000 

5003856000 

35026992000 

17513496000 

5837832000 

35026992000 

245188944000 

1915538625 

r 1 

0 

-171811 

-2566373 

3468455 

301286857 

-4590637187 

-14398288259 

421762048 

0 


22861440 

38102400 

4191264 

62868960 

1089728640 

1089728640 

127702575 

r 2 

-16459 

-9089 

-14339 

-163307 

-703993 

-2987689 

56979607 

-23168 

0 


30240 

17280 

51840 

285120 

95040 

673920 

3991680 

45045 

r 3 

0 

-64649 

-1946191 

-3764137 

32748677 

101232611 

-17097977 

-2091088 

0 

3265920 

5443200 

2993760 

8981280 

11975040 

2395008 

1403325 

r 4 

81491 

141751 

288599 

6701891 

560257 

-4187303 

78901349 

642734 

0 

453600 

777600 

777600 

4276800 

1425600 

777600 

59875200 

467775 

r 5 

0 

88517 

61633 

-93301 

-1620853 

1038857 

255833 

-538927 

0 


5443200 

3024000 

151200 

1360800 

604800 

604800 

850500 

r 6 

-11737 

-119603 

-269467 

16957 

218939 

-488741 

-5818427 

773411 

0 

340200 

2332800 

2332800 

1166400 

388800 

2332800 

16329600 

4082400 

r 7 

143 

298727 

630773 

1057927 

-1137851 

-152395 

1823393 

-1825543 

0 


145152 

45722880 

15240960 

15240960 

9144576 

3048192 

15240960 

45722880 

r 8 

3223 

28127 

-7337 

-9859 

7069 

157289 

-966457 

58621 

0 

793800 

5443200 

5443200 

388800 

907200 

5443200 

38102400 

9525600 

r 9 

-143 

-30173 

-79651 

4771 

61009 

-44227 

24421 

-69019 

0 

435456 

19595520 

32659200 

1306368 

19595520 

6531840 

6531840 

97977600 

r io 

-3211 

-403 

7033 

91 

-247 

11519 

-32227 

2447 

0 

13608000 

11664000 

11664000 

2916000 

243000 

11664000 

81648000 

40824000 

r ii 

13 

13 

-13 

-221 

13 

-221 

53 

-299 

0 


483840 

207360 

345600 

2280960 

84480 

2280960 

1774080 

79833600 

r 12 

1 

-1 

-1 

7 

-1 

1 

-19 

1 

0 


181440 

155520 

155520 

427680 

71280 

155520 

11975040 

5987520 

r 13 

-1 

-1 

29 

-13 

113 

-173 

1 

-47 

0 


1451520 

2612736 

21772800 

9580032 

143700480 

622702080 

17791488 

9340531200 

r 14 

-1 

1 

-1 

1 

-1 

1 

-1 

1 

0 


25401600 

10886400 

10886400 

17107200 

39916800 

141523200 

838252800 

10897286400 

r 15 

1 

-1 

1 

-1 

1 

-1 

1 

-1 

0 

203212800 

261273600 

435456000 

958003200 

2874009600 

12454041600 

87178291200 

1307674368000 


B. Inversion of the IBSE-A: system (18) for C — A 

For the special case of a Fourier discretization of the Poisson problem, where C is the periodic Laplacian, 
the method of inversion given in Section 3.4 is complicated by the nullspace of A. In particular, the main 
block in Equation (18), 


( A -xe A\ 

V Hk ) 


(60) 


is not invertible, and thus the Schur-complement SC (32) cannot be directly formed and prefactored. The 
resolution of this problem is conceptually simple but computationally involved, so we will give a detailed 
description for the simpler case of a direct-forcing Immersed Boundary discretization and state the result for 
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our problem. A direct-forcing IB discretization of the Poisson problem with Dirichlet boundary conditions 
is 



(61) 


Formally, G may be computed by solving 

{S*A~ 1 S)G = S*A~ 1 f - b. (62) 

Despite the fact that Equation (61) is invertible, the periodic Laplacian has a nullspace, making direct 
implementation of this approach to compute G impossible. Instead, let us decompose the solution u as 


u = uq — Uv, (63) 

where uq has mean 0, U is a scalar, and v = l(Ax) d , where 1 denotes the vector of all ones. We note that 
v spans the nullspace of the self-adjoint operator A, and we have scaled v so that v T u is equal to the the 
discrete integral of u. We will also denote the averaging operator f c w = I Cl 1 f c w, discretely this is given 
for C by f c w = (2n)~ d v T w and for T by f r H = {2 tt)~ 1 As J H, where As is the vector of quadrature weights 
on the boundary (see Section 3.1). Plugging the decomposition given in Equation (63) into Equation (61), 
we find that 


Au 0 + SG = f. 

Define a family of pseudo-inverses to A, {-4 M } M >o, by its Fourier series to be 


(64) 


T'f-jvfo 1*1 = 0 , 
^ 1 l*l*o. 


(65) 


Because uq has mean 0, A^Auq = uq for any ft. Applying A M to Equation (64), we find that uq = 
Apf — A^SG, and applying S* to this equation, along with the constraint that S*u = 6, yields that 


S*A I1 SG + S*(Uv) = S*A fJi f - b, 


( 66 ) 


which is a reformulation of the formal equation (62). This is now an equation for the two unknowns G 
and [/, and must be supplemented with the additional equation found by projecting Ait 4- SG = f onto 
the nullspace of A, giving v T SG = v T f. The combination of these two equations forms an augmented 
Schur-complement system 


[S*A tl S S*v\ (G\ _ (S*A ll f - b\ 
\ id .S' ) \U) \ v\f J 


that may be directly formed and prefactored. With (G, U) known, we may compute Uq = A^f — A^SG, 
and finally it = uq — Uv. This system may be simplified by exploiting the fact that the nullspace of A is 
only constant functions and that f c SG = f r G, giving the equivalent system 


(S*A m S {3- 1 !) (G\ _ (S*A lt f - b\ 
{AsT/2n ){UJ f c f )' 


( 68 ) 


where f3 is a free parameter. We remark that when As is a constant vector, choosing /3 = n^dy symmetrizes 
the matrix in Equation (68). Although any choice of /i and (3 may be used, we have observed empirically 
that /jl = 0 and f3 = ribdy yields the most well-conditioned system. Once G and U are known, we may 
compute u to be 


u = A ll {f-SG)-p- 1 U. 
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(69) 


Following the same procedure for Equation (18) gives the augmented Schur-complement: 


T* k A,S p-'Y h \ (F\ (S*A llX af-a\ 

S*A^XE^{'H k )~ l T k S*A fl S G = T*A llX nf . (70) 

V f G XEm k T l T k f c S ) \xj V fcXnf ) 

Here Y k is defined by Y k = (l nbdy ®(fc-i)n b d y ) T ; this f° rm is due to the fact that the estimates of values 
produced by T k are affected by changes in mean while estimates of normal derivatives are not. Again, we 
empirically observe that choosing fj, = 0 and ft = nbdy provides the most well-conditioned system. Once F 
and G are known, we can compute f as 


Z = -{U k )~ 1 T k F, 

(71) 

and once £ is known, u can be computed as 

u = A^xnf + XeA£ - SG ) - P~ X U. 

(72) 
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